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1.0 Introduction 


NASA’s Ice, Cloud, and Land Elevation Satellite (ICESat) mission will be launched late 2001. 
It’s primary instrument is the Geoscience Laser Altimeter System (GLAS) instrument. The main 
purpose of this instrument is to measure elevation changes of the Greenland and Antarctic ice 
sheets. To accurately measure the ranges it is necessary to correct for the atmospheric delay of the 
laser pulses. The atmospheric delay depends on the integral of the refractive index along the path 
that the laser pulse travels through the atmosphere. The refractive index of air at optical wave- 
lengths is a function of density and molecular composition. For ray paths near zenith and closed 
form equations for the refractivity, the atmospheric delay can be shown to be directly related to 
surface pressure and total column precipitable water vapor. For ray paths off zenith a mapping 
function relates the delay to the zenith delay. The closed form equations for refractivity recom- 
mended by the International Union of Geodesy and Geophysics (IUGG) are optimized for ground 
based geodesy techniques and in the next section we will consider whether these equations are 
suitable for satellite laser altimetry. 

To estimate surface pressure and precipitable water vapor, numerical weather models are appeal- 
ing because they are internally consistent and provide spatially uniform coverage. Values for sur- 
face pressure and precipitable water vapor will be calculated from global atmospheric analyses. 
We will use the National Center for Environmental Prediction (NCEP) global numerical weather 
analyses. As of January 2000, these analyses are produced on a 1 by 1 degree grid every 6 hours. 
Fields included are temperature, geopotential height, and relative humidity at standard upper 
atmospheric pressure levels. These atmospheric fields will be interpolated to the location and 
time tag of the laser footprints. The NCEP provides a surface pressure field, however error in this 
field make it unsuitable for our purposes. Surface pressure will be calculated by integrating the 
upper atmospheric field down to the surface height as given by the initial laser footprint location. 
Appendix A describes this procedure in detail. 

The GLAS single shot error budget assumes less than 20 mm rms error in the atmospheric delay. 
For all the steps in the atmospheric delay algorithm we require an estimate of the associated error. 
Unfortunately, NCEP does not provide formal error estimates for its numerical weather models. 
The atmospheric delay algorithm will be validated using Automatic Weather Station data in polar 
regions, in-situ meteorological data, where available, at GPS sites, and the estimated delay values 
from global GPS data. The validation studies performed so far predict less than 12 mm rms error 
in delay. 

Note that all units in this chapter are SI unless otherwise stated. 


2.0 Algorithm Description 

2.1 Group Refractivity Models 

Pulsed satellite laser altimeter systems use the round trip pulse travel time to estimate the distance 
to the Earth’s surface. Atmospheric refraction increases the propagation time and this delay may 
be corrected for by integrating the refractive index along the ray path through the atmosphere. 
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The one-way correction to the GLAS range measurement, AL, due the refractive effects of the 
Earth’s atmosphere is given by 



( 1 ) 


where n(s) is the group refractive index of the atmosphere along the ray path, S ATM is the curved 
path followed by the laser pulse from the space craft to the ground, and S VAC is the straight line 
path from the space craft to the ground. Evaluation of the second integral only requires the space 
craft and laser footprint coordinates. Evaluation of the first integral also requires knowledge of 
the refractive index along the curved ray path and is most accurately calculated using ray tracing 
or numerical integration methods. This is not practical for large amounts of data, models that 
relate the total delay to the zenith delay by a mapping function are commonly used such that 
[Davis, et al., 1985] 


where m(e, P) is a mapping function that depends on elevation angle, e , and a parameter vector 
P, and the integral is evaluated along a zenith path from ground point, Z, to the space craft to give 
the zenith delay. The mapping function will be investigated in a following section. 

The argument of the integral, (n - 1 ) is the refractivity and is normally given in parts-per-million 


sure, temperature, and composition. Precision geodetic surveying requires an accurate knowledge 
of the group refractive index of air. The first comprehensive set of equations was given by Elden 
(1953), these equations were the standard recommended by the International Union of Geodesy 
and Geophysics (IUGG) for optical wavelengths until recently [ IUGG , 1963], Other formulations 
have been used as new empirical measurements of refractivity have demonstrated errors in the 
original Elden equations [Elden, 1966; Owens, 1967; Peck and Reeder, 1972; Birch and Downs, 
1994] . A working group of the IUGG has updated their recommendations for the refractive index 
of air at optical frequencies [IUGG, 1999]. There are two solutions recommended, one is a simple 
closed formula with an accuracy of no better than one part per million and the other is a more 
complicated algorithm when better that one part per million accuracy is required. 

The more accurate algorithm is documented in Ciddor [1996] and Ciddor and Hill [1999]. The 
Ciddor formulas incorporate the latest moist air refractivity data. They match the measurements 
within experimental accuracy and are expected to be accurate over a wide range of atmospheric 
parameters to a few parts in 10 8 . The density equations used are valid over ranges of at least -40 to 
100 °C, 800 to 1200 hPa, and 0 to 100% relative humidity. However, the algorithm requires a 
number of steps and as such is unsuitable for practical integration through the atmosphere. It may 
be used as baseline measurement to compare closed form approximations of refractivity. 

The simple closed formula recommended by IUGG is the same as the formula recommended in 



( 2 ) 


i.e., (n - 1 ) = 10 6 N . Refractivity varies in the atmosphere primarily as a function of pres- 
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1963, amended for an increased C0 2 value of 375 ppm from 300 ppm. This value of carbon diox- 
ide is expected to be valid in 2004 based on estimations of rate of increase and is recommended as 
the standard value to be used in all formulas. The formula for refractivity is 


N 

N s 


"’Mi 


4.88660 

287.6155 + — + 


T 

0.06800 


(3a) 


where N is group refractivity, T is temperature (K), P is total pressure (hPa), P w is water vapor par- 
tial pressure (hPa), N s is standard air refractivity with 375 ppm C0 2 , T = 273.15 K, P = 
1013.25 hPa, P w = 0 hPa, X is wavelength (pm). For the GLAS wavelength of 1 .064 pm, and con- 
verting pressure units to Pa from hPa, equation 3a becomes 

P P 

N = (0.7871275 K/Pa)^ - (0.1127 K/Pa)^ (3b) 


The IUGG closed formula is discussed in Riieger [1996], it was noted that the humidity term is 
not very accurate and may have other systematic errors. Riieger strongly recommended that the 
equations in Owens [1967] be used for higher precision measurements. For optical frequencies, 
the group refractivity is given by 


N = k x (X) 



+ k 2 (X)^Z 


-l 

w 


k x {X) 

k 2 (X) = 


= 164.63860 < 238 - Q185 + ^ + 4.77299 (57362 + 

(238.0185 - X 2 ) (57.362 - X 2 ) 

0.648731 + 0.0174174X -2 + 3.55750x10“ 4 ?i“ 4 + 6. 1957x10“ V 6 


(4a) 


where k/T.) and k 2 (X) are experimentally determined functions of the laser wavelength (K/Pa), X 
is wavelength (pm), P d and P w are the partial pressures of dry-air and water vapor (Pa), T is tem- 
perature (K), and Z d and Z w are the compressibilities of dry-air and water vapor. The Owens equa- 
tions are based on 300 ppm C0 2 . By modifying equation (2) from Ciddor [1996], a correction 
factor for dry air refractivity based on carbon dioxide concentration may be calculated such that 


F 


c 


Nl 

*1 


1 + 


(c 2 ~ c l) 

(q + 1.8722X10 6 ) 


(5) 


where c 2 and c 2 are carbon dioxide concentrations in ppm for the respective dry air refractivity val- 
ues. For our application, c, is 300 ppm and c 2 is 375 ppm, this gives a correction factor of 
F c = 1 .000040053. Note that this is only applied to the first term in Equation 4a. For the GLAS 


5 


wavelength of 1 .064 pm, k, = 0.7866070 K/Pa and k 2 = 0.6644364 K/Pa, by also applying the cor- 
rection for carbon dioxide concentration, equation 4a becomes 

p p 

N = (0.7866385 K/Pa)yZ^ + (0.6644364 K/Pa) (4b) 

To approximately compare equation 3b to equation 4b, we will ignore the compressibility factors. 
Then equation 4b needs re-written in terms of P d = P - P w . Equation 4b then becomes approxi- 
mately N = (0.7866385 K/Pa)P / T - (0.1222021 K/Pa )P W /T . It can be seen that the two 

equations are roughly the same. Percentage wise, the constant in front of the water vapor term is 
the most different, however the water vapor partial pressure is a fraction of the total pressure. The 
following describes a more quantitative comparison of the two refractivity equations over a range 
of atmospheric conditions. 

The IUGG precision algorithm (Ciddor99) will be used as a baseline to compare the IUGG closed 
formula (IUGG99) to the modified Owens formula (Owens375). All comparisons will be carried 
out at a wavelength of 1 .064 pm, as appropriate for the GLAS laser altimeter. Setting pressure to 
1000 hPa, refractivity is compared over a wide range of temperature and relative humidity and 
shown in Figure 1 . It can be seen that the IUGG99 model performs best for temperatures above 
0 °C and high humidities. This is not surprising considering that the driving application for the 
IUGG recommendations is ground based laser geodetic surveying. The modified Owens model 
performs consistently better for lower temperatures, this suggests that Owens375 is the better 
model for integration through the atmosphere as the temperature drops through the troposphere. 


A: Owens375 - Ciddor99 B: UGG99 - Ciddor99 




Temperature (C) Temperature (C) 


Figure 1 . Comparison between closed form models of refractivity and Ciddor99 model, for a con- 
stant pressure of 1000 hPa, A, = 1 .064 pm. 
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To assess how the two closed form models perform for satellite laser altimetry we must integrate 
the differences through the atmosphere. As a point of comparison, we used standard atmospheric 
profiles for temperature and pressure [NOAA, 1976] and varied the relative humidity. The profiles 
of refractivity were numerically integrated as per Equation 2, setting the mapping function to one 
to give a zenith delay difference. Figure 2 show profiles for three different values of relative 
humidity. For zero humidity the modified Owens model is clearly preferable, it becomes less 
obvious at higher humidities. The integrated difference has been calculated over the whole range 
of humidities, Figure 3 compares the absolute value of these delay differences. Except for humid- 
ities above 80%, the modified Owens model is the preferred model. While the true atmospheric 
profile will vary from the standard atmosphere and relative humidity is almost certainly not con- 
stant along the profile, for satellite altimetry applications the preferred model appears to be 
Owens375 and will be used for the ICESat mission. Based on the comparison in Figure 3, we can 
expect inaccuracies in the refractivity model to contribute less than 1 mm to the delay error. It 
should be noted here that satellite laser ranging (SLR) stations commonly use a refractivity model 
that is the same form as IUGG99, combined with specialized mapping functions [Marini and 
Murray, 1973], this is the International Earth Rotation Service standard model. We suggest that 
the standards might be updated to account for the fact that the models currently used are best 
suited for ground based observations. 
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A: Comparison for Standard Atmosphere, RH = 0% 



C: Comparison for Standard Atmosphere, RH = 100% 



Figure 2. Comparison between closed form models of refractivity and Ciddor99 model, for stan- 
dard atmosphere profiles of temperature and pressure and varying relative humidity, 
X = 1 .064 Jim. 


0.6 



Figure 3. Absolute integrated delay differences between closed form refractivity models and 
Ciddor99 model, for standard atmosphere profiles of temperature and pressure, X = 1 .064 pm. 
Dotted vertical line marks the cross-over point for best model with respect to absolute delay dif- 
ferences. 

2.2 Zenith Delay Equations 

The modified Owens refractivity equations use compressibilities to account for non-ideal gas 
behavior such that 


Pi 


T R ' 


( 6 ) 


where p, is the density of gas i (dry air or water vapor in the case of the atmosphere) with molecu- 
lar weight M. and compressibility Z , at pressure P and temperature T; and R is the Universal gas 
constant. For reasons that will become clearly shortly, Equation 4 is often combined with Equa- 
tion 6 and written such that the total density of gas, p = p rf + p w , appears. In this form Equa- 
tion 4 becomes 


N 






MJM V 


(7) 


where we have dropped the dependence of k 1 and k 2 on wavelength for simplicity and the dry air 
refractivity correction factor for carbon dioxide is F c , as given in Equation 5. The reason for this 
choice of form for refractivity is that the first term is the largest term and with the assumption that 
the atmosphere is in hydrostatic equilibrium the integral in Equation 2 for the range correction can 
be solved exactly. The assumption of hydrostatic equilibrium is violated only under extreme 
weather conditions, such as thunderstorms and heavy turbulence, where there are significant verti- 
cal accelerations [ Fleagle and Businger, 1980], These accelerations can reach 1% of gravity, 


9 


which correspond to a delay error of approximately 20 mm [ Davis et al., 1985], however these 
extreme conditions are rare and the GLAS laser altimeter will not range through thunderstorms. 

We use the hydrostatic equation 


JZ = -p(z)g(z) (8) 

where z is height through the atmospheric column and g(z) is the gravitational acceleration. 
Substituting Equation 8 into the first term of Equation 7 and then substituting into Equation 2 
yields the “hydrostatic” component of the zenith range correction, AL^, which can be written as 

oo 

AL h = (9) 

d z 

where g m is the mean value of gravity in the column of the atmosphere. Since gravity decreases 
slowly with height and can be closely approximated as a simple function of latitude, this value can 
be expressed accurately in terms of the height, Z, and latitude, ()) , of the ground point to which the 
altimeter measurement is made [Saastamoinen, 1972] 


g m = 9.8062(1 - 0.00265 cos (2^) - 3.1x10 7 (0.9Z + 7300)) ms 2 (10) 

Equation 9 can be further reduced because the integral is simply the surface pressure at height Z, 
and therefore the largest part of the atmospheric range correction in the zenith direction is given 
by 


a l h = ^F c k± g -y s 


( 11 ) 


where P s is the surface pressure. 

The remaining part of the zenith atmospheric range correction is due to the residual part of the 
water vapor not included in the hydrostatic term, commonly called the “wet” component, A L w . 
Substituting the second term of Equation 7 into Equation 2 gives 


M w 

where k 2 = k 2 - F c k ] . The integral is simply the total column precipitable water vapor, 

PW, an atmospheric variable often reported in atmospheric models. The zenith wet delay can now 
be written as 


A L w = 10 -%'^-PW (13) 

When the empirical functions in the refractivity equation given in Owens [1967] are evaluated for 
the GLAS laser altimeter operational wavelength of 1.064 |am, k ( = 0.7866070K/Pa and 

k 2 = 0. 6644364 K/Pa . Dry air with a carbon dioxide concentration of 375 ppm has a molecular 

weight of M d = 28.9632 kg.kmol 1 [Ciddor, 1996], and water has a molecular weight of 
M w = 18.0152 kg.kmol ‘. Using these values and the previously calculated value for the carbon 
dioxide correction factor of 1.000040053, k 2 = 0.175 1 448 K/Pa . Combining these values and 

using the Universal gas constant value of R = 8314.510 J-kmoU.K 1 into Equations 11 and 13 
gives the final zenith delay equations 

AL Z = AL h + AL w (14) 

A L h = (2.2582rn 2 s 2 /Pa)g m 1 P 5 (15) 

A L w = (8.0834xl0“ 5 m/mm)PW (16) 

Given an average surface pressure value of 1000 hPa and an approximate value of 9.8 ms 2 for the 
mean gravity, the zenith hydrostatic delay is approximately 2.3 m and is the major component of 
total delay. Zenith wet delay is much more variable, given precipitable water vapor values of less 
than 10 mm in the polar regions to 50 mm in the tropics, it varies from 1 to 4 mm. 

2.3 Off-Nadir Pointing Corrections 

The mapping function relates the total atmospheric delay at an arbitrary elevation angle to the 
zenith delay such that 


A L = m(e, P)AL Z (17) 

where e is the elevation angle and P is a vector that commonly consists of various climatological 
parameters. The mapping function assumes azimuthal symmetry of the atmosphere about the 
ground point, a very good assumption for the near-nadir pointing ICESat mission, although hori- 
zontal gradients are a significant error source for low elevation angle satellite laser ranging. 

When it is assumed that the refractivity of the troposphere is azimuthally and spherically symmet- 
ric, Marini [1972] showed that the continued fraction form of the mapping function is 
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m(e) = 


1 


a 


(18) 


sine + 

b 

sine h 

c 

sine + — 

sine + ... 

where a, b, c, etc., are parameters that may be approximated using climatic data. The very sim- 
plest form of this equation is 


m(e) 


1 

sine 


(19) 


which works best at elevation angles near zenith since that parameters a, b, c, etc., are all signifi- 
cantly less than one. A number of different forms of the mapping function have been published, 
to test the accuracy of using Equation 19 we will compare this to two different but widely using 
mapping functions. One is by Davis, et al. [1985], named CfA-2.2, which depends on surface 
pressure and surface temperature. The other is by Niell [1996] which depends on latitude and day 
of year. The different climatic variables used are due to the different climatologies and functional 
forms of the parameters used. We compared the simple mapping function to the test functions by 
subtracting the test function from the simple mapping function and multiplying by 2.3 m, which is 
a rough estimate for zenith delay. This gives estimates of how much the total delay will change 
by. Figures 4 and 5 show these comparisons over a range of climatological conditions. 
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Figure 4. Change in delay of the simple mapping function compared to CfA-2.2 mapping func- 
tion. Left plot is for P 0 = 1000 mbar, right plot is for T 0 = 0°C. 
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Change in Delay (mm) 



Change in Delay (mm) 



Figure 5. Change in delay of the simple function compared to Niell mapping function. Left plot 
is for maximum day of year phase, right plot is for minimum day of year phase. 


For both of the comparisons, using the simple form of the mapping functions compares very 
closely to the other forms. We don’t expect that the GLAS space craft to point beyond 10° off 
nadir, the differences in this region are less than 0.5 mm for CfA-2.2 and 0.1 mm for Niell. It 
should be noted that these other mapping functions are optimized for low elevation angles and in 
fact we expect the Niell mapping function to work better at higher elevation angles due to its func- 
tional form. So we will use the simple m(e) = l/(sine) form of the mapping function. The 
error associated with the mapping function should be less than 0.5 mm for pointing angles of less 
than 10° and decreases to zero for nadir pointing. 

There is another concern for off-nadir pointing of the space craft, which is the change in expected 
footprint location due to bending of the ray in the atmosphere. This effect will not significantly 
change the atmospheric delay calculation but should be considered for spacecraft pointing cali- 
brations where the location of the laser footprint is directly measured at the ground. The real 
curved path is shown by the dashed line in Figure 6. Pj is the expected ground location of the 

laser footprint for the satellite position and pointing angle, a, = 90° - £ ( , as measured at the 

satellite. P 2 is the real ground location of the laser footprint after following the refracted path 
through the atmosphere, which is shifted by a distance d towards the sub-satellite point. If the sat- 
ellite position and real footprint location were used to calculate the apparent satellite pointing 

angle, a 2 = 90° - e 2 , this would be in error by a certain amount such that a, = a 2 - 5a . 

This correction can be approximated by a simple expression for pointing angles of less than 75° 
[Astronomical Almanac , 1999] such that 
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8a = 0 .00452 ° P tan a 2 /( 273 + T) 


( 20 ) 


where T is the temperature (°C) and P is the pressure (hPa) at the surface. Using standard surface 
values of 15°C and 1013 hPa, this gives an approximate value of 


8a = 0.016° tan a 2 
= 57 tan a 2 


( 21 ) 


At an altitude of 600 km and pointing angle of 10 degrees, the pointing error will be approxi- 
mately 10 arcseconds and the distance the laser footprint is shifted by will be 30 m. 



Figure 6. Geometry of laser ray path, dashed line is real path through atmosphere. 


As an aside, this correction equation may also be used to estimate whether horizontal gradients in 
the pressure fields will greatly affect the path of the laser pulse. A typical upper value of the syn- 
optic pressure gradient is 10 mbar per 100 km. Near the surface the derivative of pressure with 
respect to height is approximately 0.1 hPa/m. This means that a typical slope to the pressure field 
is 2 arcseconds. Assuming that this gradient is constant through the atmosphere (actually, should 
decrease exponentially) will can put this angle into Equation 21 to see how much the path will 

-4 

deviate. For a 2 arcsecond slope the deviation is 6x10 arcseconds, which is entirely negligible. 

2.4 NCEP Global Analyses 

The zenith delay formulas given in a previous section are directly dependent on surface pressure 
and total precipitable water vapor. We require a data set that will allow us to calculate values for 
surface pressure and total precipitable water vapor at the laser footprint locations. Numerical 
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weather models are appealing because they are internally consistent and provide spatially uniform 
coverage. We will use the global analyses and forecasts produced by the National Oceanic and 
Atmospheric Administration (NOAA) National Center for Environmental Prediction (NCEP), the 
data products are from the Global Data Assimilation System (GDAS) [Kanamitsu, 1989; Kana- 
mitsu et. al., 1991]. 

The NCEP GDAS uses a spectral model based on the primitive atmospheric equations, observa- 
tional data is assimilated using a spectral statistical-interpolation (SSI) analysis [ Parrish and Der- 
ber, 1992], Data sources used to create the NCEP global analyses include ground stations, 
radiosondes, satellites, and buoys. They are produced on 1 degree uniform latitude and longitude 
grid every 6 hours, starting at 0 GMT. These analyses and forecasts consist of a number of mete- 
orological fields for a standard set of levels from the surface to the stratosphere. The fields that 
we will use for our surface pressure model are temperature, geopotential height, and relative 
humidity for the tropospheric pressure levels between 1000 mbar and 300 mbar. The total precip- 
itable water vapor is given as a single field integrated through the entire atmospheric column and 
may also be calculated by integrating relative humidity through the tropospheric levels. 

The NCEP analyses report a surface pressure field, however it is unsuitable for our purposes. 
The topography field it is produced on has large errors and the spectral interpolation used creates 
spurious waves in the field near areas of rapid change in elevation, as encountered at the edges of 
ice sheets. It is essential to perform this integration rather than use the NCEP surface pressure 
field so that the correct height of the laser footprint is used. We have compared the NCEP surface 
pressure field with ground station data in Antarctica and found it to be biased by as much as 40 
mbar. An atmospheric model of pressure with respect to height is required to reduce the upper 
level NCEP fields to a surface pressure value. A hydrostatic equilibrium model of the atmosphere 
is integrated from an upper atmospheric level to the estimated height of the laser footprint in order 
to calculate surface pressure. This process is described in more detail in the next section. The 
NCEP global analyses give total column precipitable water vapor as a single field evaluated at the 
surface, we will use this without modification as input into the zenith wet delay equation. The 
precipitable water vapor contribution to total delay is small but highly variable both spatially and 
temporally and should be monitored throughout the ICES AT mission. 

The NCEP global analyses can be downloaded near real time from an anonymous NOAA ftp site 
(ftp.ncep.noaa.gov). It should be noted that this ftp site is not an archive and the products only 
remain available for approximately 24 hours. Throughout the mission the ICESat team will 
download and archive the NCEP fields that is required to calculate atmospheric delay. 

There are two different runs of the GDAS we will use, based on availability: (a) the final run 
(FNL); and (b) the aviation run (AYN). The final run produces the best analysis, as it is delayed to 
allow for late arriving data. The analysis and a 6 hour forecast are posted to the ftp site approxi- 
mately 6 to 10 hours after the analysis time, the 00Z and 12Z analyses take longer than the 06Z 
and 18Z analyses because there are more ground station data taken on 12 hour intervals. The avi- 
ation run uses exactly the same model as the final analysis except that it is run at an earlier time 
and therefore has less data included, the AVN analysis is posted to the ftp site approximately 3.5 
hours after the analysis time. The forecasts are every 6 hours out to 84 hours and take an addi- 
tional 5 minutes to be posted for each forecast, all forecasts are posted to the ftp site approxi- 
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mately 5.5 hours after the analysis time. Table 1 is a summary of approximate NCEP delivery 
times. 

Over 95% of the time, the final analysis is available for use. However, sometimes this run is not 
performed by NCEP due to time constraints. Our next best choice is the aviation analysis, which 
is almost always available. If the analyses are not available we will use the forecasts, the order of 
preference being the FNL 6 hour forecast, then the AVN forecasts in ascending time. The atmo- 
spheric delay estimates should be flagged according to the data used, in order to assess the quality 
of the estimate. The forecasts need to be archived until it is certain we can obtain the better anal- 
yses. 


Table 1: NCEP Delivery Times as of March 29, 2000 


Model 

OOOOZ 

0600Z 

1200Z 

1800Z 

FNL 





analysis 

0730Z 

1030Z 

2130Z 

2300Z 

6 hr forecast 

0735Z 

1035Z 

2135Z 

2305Z 

AVN 





analysis 

0315Z 

0920Z 

1515Z 

2115Z 

6 hr forecast 

0325Z 

0930Z 

1525Z 

2125Z 

84 hr forecast 

0405Z 

1010Z 

1605Z 

2205Z 


The NCEP products are stored in GRIB format (GRIdded Binary). This format is widely used in 
the meteorological community and is the World Meteorological Organization standard for 
exchanging gridded binary data. It is thoroughly described by NCEP Office Note 388 [ Stackpole , 
1994] . NCEP has codes for reading GRIB format that we have incorporated into our atmospheric 
delay software package. 

Unfortunately, there are no formal error estimates provided for the atmospheric fields produced by 
the NCEP GDAS. Studies have compared analyses produced by different forecasting centers, 
notably the European Center for Medium-Range Weather Forecasting (ECMWF), however these 
competing analyses use much the same input data and physical models and as such do not provide 
a quantitative error estimate [Trenbreth & Olsen, 1988; Boer et al., 1992] . Validation studies have 
been performed for the surface pressure and precipitable water vapor estimates to address this 
shortfall. 

2.5 Surface Pressure 

An atmospheric model of pressure with respect to height is required to reduce the upper level 
NCEP fields to a surface pressure. To simplify the physical model of the atmosphere we will 
make certain assumptions. A static atmosphere model will allow us to consider the vertical distri- 
bution of atmospheric variables. Although the atmosphere is actually a dynamic system, static 
atmosphere formulas for variables like pressure and density are valid to a high degree of accuracy. 
We will assume a horizontally stratified atmosphere in hydrostatic equilibrium, such that pressure 
is related to height by the hydrostatic equation 
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dP = -g(Z)p(Z)dZ 


( 22 ) 


where Z is geometric height, P is pressure, g is gravity, and p is density. 

To allow easier integration of this equation, we will convert geometric height into geopotential 
height. A geopotential meter is defined as the work done by lifting a unit mass one geometric 
meter through a region in which gravity is uniformly 9.80665 m/s 2 , the value of mean sea level 
gravity. The geopotential measured with respect to mean sea level (assumed zero potential) is 
called geopotential height, H, such that 


H = - fgdz. (23) 

#o J o 

where g 0 = 9.80665 m 2 /s 2 m’ [NOAA, 1976], The derivative of this equation with respect to geo- 
metric height is 


g 0 dH = gdZ (24) 

This can be substituted into the hydrostatic equation to give 

dP = -g 0 p(H)dH (25) 

We now require an expression that will convert elevation in geometric meters to geopotential 
meters. This will be related to the variation of gravity with height. Approximating the Earth as a 
sphere with only radial mass variations, gravity is inversely proportional to radius squared, which 
will give a conversion equation of 


8 


8 msl 


R 


(R + zy 


(26) 


where R = 6371009 m is mean radius of the Earth, g msl is gravity at mean sea level. Substituting 
this equation for gravity into Equation (23) gives the conversion formula 


H 


8 msl 


RZ 


g o (R + z ) 


(27) 


Mean sea level gravity depends on geodetic latitude, the formula is based on calculations of the 
standard geodetic refence system [Moritz, 1980] such that 


_i 

8msl = g e(? (l + £sin 2 <|))(l - e 2 sin 2 <|)) 2 (28) 
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where <|> is latitude, g eq = 9.7803267715 m/s 2 , k = 0.001931851353, e 2 = 0.00669438002290. 


Many atmospheric models, such as the U.S. Standard Atmosphere [NOAA, 1976], simplify their 
calculations for pressure by assuming the air to be a dry, ideal gas. We shall include non-ideal gas 
effects and water vapor partial pressure. The equation of state for a pure non-ideal gas is 


„-i PV _ m 

RT ~ M 


(29) 


where Z 1 is called the inverse compressibility and depends empirically on pressure and tempera- 
ture [Harrison, 1965b], P is pressure, V is volume, R is the universal gas constant, T is tempera- 
ture, m is mass, and M is molecular weight. 

Density can be written as p = m / V , and we can split mass components of water and dry air: 
m = m w + m d . If we assume that moist air obeys Dalton’s Law of partial pressures, the sepa- 
rate masses can be evaluated by the non-ideal equation of state to give a density equation of 

p = ■±f(Z-JP H ,M w + Z- J i (P-P w )M d ) (30) 


where R = 8314.510 J/kmol.K, M w = 18.0152 kg/kmol, M d = 28.9632 kg/kmol for 375 ppm carbon 
dioxide concentration, P is the total pressure and P w is the partial pressure of the water vapor in 
the air. It is implicitly assumed that the dry air components are homogeneously mixed throughout 
the lower atmosphere and therefore the mean molecular weight of dry air is a constant. 

Equations for inverse compressibility have been experimentally determined by Owens [1967], 
these formulas are accurate to within a few parts per million 

p 

Z~J = 1 + 1650— ^[1 - 0.01317(7 - 273.15) + 1.75xl0“ 4 (r - 273. 15) 2 (31) 

T ^ 

+ 1.44xl0“ 6 (r - 273.15 ) 3 ] 



1 +(P~P W ) 


57.90x10 1 + 


0.52 


9.4611x10 


,-4 (T - 273.15) ' 
T 2 


(32) 


We need an equation for water vapor pressure. Given sufficient saturation vapor pressure data 
over a wide range of temperature, the information can be stored in an analytical form. One of the 
better forms uses Chebyshev polynomials [McGarry, 1983] 
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( 33 ) 


P ^ 

7 ' l °giofp £ ) = T + I "A<*> 

v h ' , . i 

where P b = 1000 Pa, P is saturation vapor pressure, E s (x) are Chebyshev polynomials: 


E 0 (x) 
E i(x) 

E s+ iW 


1 


X 

2 xE s (x) - E s _ x (x) > 


2T - (T + T ) 

v max mm' 

T - T • 

max mm 


(34) 


The coefficients a s (5 = 0,...,10) are a = {27 9 4 .027 , 1430.604, -18.234, 7.674, -0.022, 0.263, 
0.146, 0.055, 0.033, 0.015, 0.013}, T max = 648 K, and T mm = 273 K [Ambrose, 1987], 

To get the water vapor pressure from the saturation vapor pressure we use relative humidity [Har- 
rison, 1965a] 


P w = RhP s (35) 

where relative humidity is in a fractional form with values between 0 and 1 . Actually, this equa- 
tion is only true for pure water vapor, not moist air. However the equation is approximately true 
for moist air. Note that the World Meteorological Organization (WMO) has adopted the practice 
of evaluating relative humidity with respect to liquid water at all temperatures, even those below 
0°C. 

We now have an expression for density that depends on temperature, relative humidity and pres- 
sure. To solve the hydrostatic equation we must express temperature and relative humidity as 
functions of geopotential height, in order to get an expression for density that only depends on 
geopotential height. The NCEP global analyses have values for temperature, geopotential height 
and relative humidity at standard pressure levels. We shall assume that temperature varies linearly 
with respect to geopotential height between these levels, a relatively good assumption for the 
lower atmosphere such that 


T = T 0 + L(H - H 0 ) 


(36) 


L 


T_ 1 Z T, 
H i ~ H o 


(37) 


where L is the temperature gradient; T 0 and H 0 are temperature and geopotential height at the 
upper level; 7} and H t are temperature and geopotential height at the lower level: H x < H < H 0 . 
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We will also assume that relative humidity varies linearly with respect to geopotential height 
between levels such that 


Rh = Rh 0 + S(H - H 0 ) 


(38) 


Rh l - Rh 0 


(39) 


where S is the relative humidity gradient, Rh 0 is relative humidity at the upper level, Rh t is relative 
humidity at the lower level. Given these expressions for temperature and relative humidity, the 
hydrostatic equation becomes 

% = ~[Z- v l (H)P v (H)M v + Z~ l (H,P)[P-P v (H)]M a ] (40) 


This differential equation is first order, non-linear and inhomogeneous, we are not able find an 
analytic solution. To obtain a numerical solution for pressure we will numerically integrate down 
from the upper level geopotential height to the desired geopotential height. Pressure varies 
smoothly with geopotential height, this means that we are relatively unrestricted in our choice of 
numerical method. We will use the Bulirsch-Stoer method, this method is one of the best ways to 
obtain high accuracy solutions with minimal computational effort, so long as integrated function 
is smooth and has no singular points within the range of integration (Press et al., 1989). 

2.6 Precipitable Water Vapor 

The NCEP global analyses give total column precipitable water vapor as a single field evaluated at 
the surface, we will use this without modification as input into the zenith wet delay equation. The 
precipitable water vapor contribution to total delay is small but highly variable both spatially and 
temporally and should be monitored throughout the ICES AT mission. 

To validate the precipitable water vapor fields we will compare them to ground station data. One 
of resources we will use is the GPS global network. Precipitable water vapor can be derived from 
estimates of GPS tropospheric delay made at each global station. This derivation requires the 
knowledge of surface pressure at the station. There are over 35 GPS stations that report surface 
pressure from on site met packages, however only two of these are in the polar regions, both in the 
north. Where directly measured surface pressure is unavailable we can use our own surface pres- 
sure model without a significant loss of accuracy. There are currently 4 stations in Antarctica and 
6 stations in the Arctic where we can make precipitable water vapor measurements. The major 
advantage of using the GPS global network is the rapid availability of data and confidence that the 
data will be available over the length of the ICES AT mission. 

2.7 Delay Correction with Respect to Height 

As can be seen in the previous sections, calculation of the surface pressure and therefore delay 
requires a knowledge of the height of the laser footprint location. The atmospheric delay correc- 
tion will be estimated early in the GLAS processing and there may be later adjustments to the 
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spacecraft orbit and footprint location height. We need a simple correction function that would be 
accurate for height changes in the range of ±100 m. 

Given that pressure varies approximately exponentially with height, we expect the correction to 
be of the form 


P' = Pe 


(-A(H'-H)) 


(41) 


where P’ and H’ are the corrected pressure and geopotential height, P and H are the original pres- 
sure and geopotential height, and A is the correction parameter. 

Neglecting water vapor, the hydrostatic equation becomes 


dP 

dH 


-8pZd PJK 

RT 


(42) 


Assuming that temperature and inverse compressibility are a constant with respect to height lets 
us solve Equation 42 for the correction factor such that 


A 


So 2 d 'M d 

RT 


(43) 


where P, T, and Zf are calculated at the original height. For typical surface values for temperature 
of 273.15 K, given g 0 = 9.80665 m 2 /s 2 .m’ [NOAA, 1976], and neglecting the compressibility fac- 
tor, A = 1.25 x 10 4 m 1 . This correction factor is defined for geopotential meters, it may be con- 
verted to geometric meters using mean sea level gravity instead of g 0 . Since atmospheric delay is 
directly proportional to surface pressure, again neglecting water vapor, this same correction factor 
may be applied such that 


A L' 


A Le 


(-A(H'-H)) 


(44) 


2.8 Spatial Interpolation 

The NCEP global analyses we will use are given on a 1 by 1 degree uniform latitude and longi- 
tude grid. We require an interpolation scheme that will allow us calculate the atmospheric field 
values at the laser footprint locations. This interpolation method will have to be computationally 
efficient to keep up with the real time data processing requirements. The global analyses have the 
highest realistic spatial resolution by design, therefore a complicated interpolation scheme 
intended for sparse data sets would not be useful nor appropriate. The upper level fields of tem- 
perature, geopotential height and relative humidity are quite smooth, a bilinear interpolation of the 
grid will be sufficient. The precipitable water vapor field is much more variable, however its 
accuracy and small contribution to the total delay do not warrant anything more complicated than 
bilinear interpolation as well. 
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Bilinear interpolation has the form 


/(4>, X) = a + bX + cY + dXY 
a = f($ v X,) 

b = /(<>!, X 2 )-/(< f»i, Xj) 

C = f(<\> 2 , X, ) - /(<!>!, Xj) 

j = /((|) l5 x^ + f(§ 2 , x 2 ) - f($i, x 2 ) - y ( <f> 2 , Xj) 

X = (X - X 1 )/(X 2 - Xj) 

F = (<) - < j ) 1 )/((|> 2 - ()>!) 


(45) 


where /is the field value, (|) is latitude, X is longitude. The subscripts 1 and 2 stand for south and 
north latitudes and the west and east longitudes of the four known grid points directly surrounding 
the unknown point. 

Any errors that arise from using bilinear interpolation will be included in the errors estimated in 
our validation studies. To estimate errors that arise purely from the bilinear interpolation, we 
decimated the 1 by 1 degree NCEP grid to 2 by 2 degrees, then interpolate back to the original 
grid. The interpolated values were then differenced from the original 1 by 1 degree grid values. 
This will give an upper bound on the interpolation errors, as we are using a coarser grid to interpo- 
late. The field we will use for testing is the 1000 hPa geopotential height field (GPH). This is the 
most appropriate field to use, considering the surface pressure algorithm integrates from an inter- 
polated upper atmospheric GPH fields down to the surface. Average global surface pressure is 
approximately 1013 hPa, therefore the 1000 hPa GPH field is the one most likely to be used in the 
surface pressure algorithm. The error in the GPH field can then be approximately converted to 
delay error using the 0.29 mm of delay to 1 meter height change correspondence. Tests indicate 
that the spatial interpolation error for the 1000 hPa GPH field, zonally averaged, is no more than 
7m, corresponding to 2 mm of atmospheric delay. 

2.9 Temporal Interpolation 

The NCEP fields used to model atmospheric delay are only produced every 6 hours. To estimate 
delay values at the times of the laser pulses we will need to temporally interpolate between the 
NCEP output times. Surface pressure is the major contributor to atmospheric delay, we can look 
at the temporal behavior of surface pressure to guide our temporal interpolation scheme. Spectral 
plots of surface pressure are used to characterize the statistical properties of the time series. For 
example, if a log-log plot of power spectral density versus frequency is has a slope of -2 then the 
time series can be described as a random walk process. For these processes, the maximum likeli- 
hood interpolator is simply a linear interpolation between adjacent points. 

Surface pressure time series from automatic weather stations (AWS) in Antarctica were used to 
calculate power spectral density plots. Stations were chosen that had 2 years of largely uninter- 
rupted time series. The AWS data used had a sampling period of 10 minutes, short gaps in the 
data were bridged was linear interpolation. Analysis of the power spectra show that log-log plots 
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of these spectra fall off at high frequencies with an approximate slope of -2.5, if some of the white 
noise tail is ignored. This slope is close to -2, therefore consistent with a random walk stochastic 
process. 

The power spectral density plots for the higher temporal resolution AWS data may be used to esti- 
mate how much variance in surface pressure we are missing by using the coarser resolution NCEP 
fields. The additional variance is the power integrated under the power spectral density plots for 
frequencies higher than the NCEP nyquist frequency. The current NCEP global analyses are pro- 
duced every 6 hours, therefore the nyquist frequency is 2 cyc/day. The additional integrated 
power was no more than 0.1 hPa for the AWS stations examined. This corresponds a delay error 
of no more than 0.3 mm due to the temporal resolution of the NCEP fields, which is negliable 
compared to the 20 mm error budget and other error sources. 

2.10 Coordinate Systems 

The calculation of atmospheric delay at the laser footprint locations requires knowledge of the 
footprint latitude, longitude, and initial height estimate. The laser footprint coordinates are given 
as geodetic latitude, longitude, and height above ellipsoid [Schutz, 1999]. The reference ellipsoid 
used is WGS-84, defined by semi-major axis (a = 6378137.0 m), semi-minor axis (b = 
6356752.3142 m), first eccentricity squared ( e 2 = 0.00669437999013) [ NIMA , 1997]. However 
all of the atmospheric data from NCEP are referenced to orthometric height, i.e. height above 
geoid. We need to convert the laser footprint ellipsoid height to height above geoid. The height 
difference can range from -100 m to 100 m, leading to a delay difference of approximately 23 mm 
to -23 mm. To calculate orthometric height we need to know the value of the geoid at the foot- 
print location. Let H be the orthometric height, h be the ellipsoid height, and N be the geoid 
value. Then the height conversion is 


H = h - N p (46) 

O 

The geoid values are calculated using the National Imagery and Mapping Agency 30 arc-minute 
geoid grid as defined for WGS-84 [NIMA, 1997], which is interpolated using bilinear interpola- 
tion to the footprint locations. 

2.11 Processing Flow 

There are two separate processes for calculating the atmospheric delays. The “background” pro- 
cess is the archiving of the required NCEP data. This data should be archived on-site as it is only 
posted to the NOAA ftp site for 24 hours before being replaced by the next day’s data. Once the 
required fields are extracted, the daily volume of data, as of March 2000, is 25.8 Mbytes. The 
main process is the calculation of the delays. The timing of the delay calculation is driven by two 
data streams: the NCEP fields and the laser ranges. We wish to calculate the delays within 12 
hours of real time, if the laser ranges are available in a timely fashion then the biggest impediment 
to delay processing is the up to 9.5 hour time lag on the 12Z FNL analysis field (Table 1). The 
following processing flow assumes that the laser range data is available on a much shorter time lag 
than the NCEP fields. All the processing is done in 6 hour time blocks, as this is the time step of 
the NCEP fields. 
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Let tj be a given 6 hourly time step, let t 2 be the next time step. Let P 12 be a set of all laser ranges 
with a time tag, t, such that t t < t < t 2 . 

1 . Between t t + 5 hr and t, + 9.5 hr: Access NCEP fields for t r If FNL analysis is archived on- 
site, proceed. If not archived, attempt to retrieve from the NOAA ftp site. If not available, use 
the next best analysis for forecast. 

2. Next: Calculate delays for P 01 at time value t r 

3. Next: Interpolate with respect to time delays for P 01 at t 0 and P 01 at t, . (Note, delays for P 01 at t 0 
were calculated in the previous loop of this process). Save these delay values. 

4. Next: Calculate delays for P 12 at time value t, . Do these calculations as soon as laser ranges 
come in, make sure completed before next time step. 

5 . Repeat for next time step (t 2 ). 

3.0 Validation 

The error budget set out in the GLAS science requirements assumes less than 20 mm single shot 
total atmospheric delay error. Total zenith delay basically consists of two components, zenith 
hydrostatic delay (ZHD) and zenith non-hydrostatic or “wet” delay (ZWD). Examples of total 
zenith atmospheric delay in the polar regions, as calculated by our model on the NCEP grid, are 
given in Figures 7 and 8. ZHD is the major component, there is approximately 2.3 mm of delay 
per millibar with a typical value of 2.3 m. Wet delay is much smaller but more variable, there is 
approximately 0.08 mm of delay per kg/m 2 of precipitable water vapor. Typical values in the 
polar regions are less than 1mm of delay and up to 6 mm of delay in the tropics. 
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Figure 8. Atmospheric delay for Greenland, 00 GMT January 1, 1998. 


To validate the atmospheric delay estimate requires that we validate our surface pressure and pre- 
cipitable water vapor estimates, especially in the polar regions. Sources of readily obtainable sur- 
face pressure measurements are Automatic Weather Stations (AWS) in Antarctica and Greenland 
[Sterns and Wendler, 1988; Steffen et al, 1996], and GPS stations with meteorological packages. 
Comparisons of model estimates with polar AWS’s for 2 year’s worth of data in 1998 and 1999 
show an rms error of less than 5 hPa once a mean offset is removed, the majority have less than 3 
hPa rms error. The offset removal is justified as the station heights are not accurately known. 
These rms surface pressure errors correspond to delay errors of less than 12 mm, usually less than 
7 mm, which is less than the 20 mm error budget. Two examples of these comparisons are given 
in Figure 9. Errors estimates are also a combination of model errors and ground station errors, the 
AWS barometers are theoretically calibrated to +/- 0.2 hPa. 
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SCHWERDTFEGER: Diff. mean = 0.4657, sdev = 1.5444, N = 2132 
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Figure 9. Schwerdtfeger AWS, Ross Ice Shelf, and Theresa AWS, West Antarctica. Upper plot 
for each station shows NCEP surface pressure values (squares) and AWS values (triangles). 
Lower plot is surface pressure differences, NCEP minus AWS. 


Precipitable water vapor data are harder to come by, however the technique of GPS meteorology 
can provide some data [Quinn and Herring , 1999]. GPS processing produces an estimate of total 
delay at microwave wavelengths, which is much more sensitive to water vapor than optical wave- 
lengths. If the surface pressure at the GPS station is known then precipitable water vapor esti- 
mates can be derived. If the station does not have an on-site met package then our NCEP surface 
pressure model may be used, greatly increasing the amount of available data. The error from 
using NCEP surface pressure is approximately 2 mm of precipitable water vapor (PW). We have 
used this method for 4 GPS stations in Antarctica and 2 in Greenland that are part of the IGS glo- 
bal network. Figure 10 shows the results of this technique, comparisons of GPS derived estimates 
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of PW versus NCEP global atmospheric model estimates of PW for a 2 week period in November 
1998. From top to bottom in the figure, the four Antarctic stations are Casey (CAS1), Davis 
(DAY1), McMurdo (MCM4), and O’ Higgins (OHIG), the two Greenland stations are Kangerlus- 
suaq (KELY), and Thule (THU1). 

First of all, it can be seen that for all the stations, except for McMurdo, the NCEP PW value is on 
average less than the GPS PW value. This may be because radiosonde humidity sensors tend to 
underestimate the true value in cases of very cold temperature and low humidities, both being true 
in the polar regions. Within the 2mm error bars for the GPS values, some NCEP values track rea- 
sonably well, i.e. Casey and Thule. The worst fitting station is O’ Higgins where the NCEP values 
are too low and too smooth. In fact for all the stations the NCEP values appear to be too smooth. 
The rms errors are all less than 4 mm of precipitable water vapor, which corresponds to 0.3 mm of 
delay, which is almost negliable and of the same order as the measurements themselves. How- 
ever, one concern is that there may be seasonal variations in the bias and error that would show up 
as seasonal variations in ice thickness. Another is that changes in the NCEP Global Data Assimi- 
lation System may show steps in the values as resolutions and data sources change. Therefore we 
will continue to model the wet delay and moniter its errors. 

The surface pressure and precipitable water vapor validation will be continuously performed and 
monitered over the lifetime of the GLAS mission. In summary, results so far indicate that the sur- 
face pressure error in polar regions is less than 5 mbar or 12 mm delay. Precipitable water vapor 
errors in polar regions are less than 5 kg/m2 or 0.4 mm delay. These errors are less than the 20 
mm error budget assumed. 
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Figure 10. Precipitable water vapor values for various stations. Solid lines with trianges are GPS 
derived values, dashed lines with circles are NCEP values. 
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Appendix A - ICESat Mission Outline 

The Ice, Cloud and land Elevation Satellite (ICESat) was launched on 13 January 
2003. The Geoscience Laser Altimeter System (GLAS) instrument onboard ICESat made 
its first laser elevation measurement of the Earth on 21 February 2003 and its last on 1 1 
October 2009. The three lasers employed by GLAS did not perform as long as expected, 
and following the failure of Laser 1 on 5 March 2003 the ICESat mission was modified to 
meet the requirement for capturing a multi-year time series of ice sheet elevations 
(Schutz et al., 2005). For the modified mission scenario, the spacecraft entered a 91 -day 
repeat science orbit (compared to a planned 183 -day repeat) and the lasers were activated 
for about 33 days of this 91 -day repeat, two or three times per year. This campaign mode 
operation is summarized in Table A.l, and other significant parameters and events are 
listed in Table A.2. ICESat laser campaigns are designated by a laser number (LI, L2 or 
L3), followed by a letter in the sequence of operation. Following campaign L2f, attempts 
to restart any of the lasers were not successful. The spacecraft was put through a series of 
engineering tests in early 2010. De-orbit maneuvers were carried out in June and July 
2010. The spacecraft was “passivated” on 14 August and reentered the Earth’s 
atmosphere on 30 August 2010 over the Barents Sea northeast of Norway. 


Table A.l: ICESat Laser Operation Campaigns 


Campaign 

Year 

Day of year 

Calendar Dates 

Number of 
days (d) 

Repeat 
orbit (d) 

Repeat tracks 1 

Lla 

2003 

051-088 

20 Feb-21 Mar 

37 

8 

001-072 to 
006-023 

L2a 

2003 

268-277/ 

277-322 

25 Sep-4 Oct/ 
4 Oct-21 Nov 

54 

8/ 

91 

028-088 to 029-100/ 
1098 to 0421 

L2b 

2004 

048-081 

17 Feb-21 Mar 

33 

91 

1284 to 0421 

L2c 

2004 

139-173 

1 8 May-2 1 Jun 

34 

91 

1283 to 0434 

L3a 

2004 

277-313 

3 Oct-8 Nov 

37 

91 

1273 to 0452 

L3b 

2005 

048-083 

17 Feb-24 Mar 

35 

91 

1258 to 0426 

L3c 

2005 

140-174 

20 May-23 Jun 

34 

91 

1275 to 0421 

L3d 

2005 

294-328 

21 Oct-24 Nov 

34 

91 

1282 to 0421 

L3e 

2006 

053-087 

22 Feb-28 Mar 

34 

91 

1283 to 0424 

L3f 

2006 

144-177 

24 May-26 Jun 

33 

91 

1283 to 0421 

L3g 

2006 

298-331 

25 Oct-27 Nov 

33 

91 

1283 to 0423 

L3h 

2007 

071-104 

12 Mar- 14 Apr 

33 

91 

1279 to 0426 

L3i 

2007 

275-309 

2 Oct-5 Nov 

34 

91 

1280 to 0421 

L3j 

2008 

048-081 

17 Feb-21 Mar 

33 

91 

1282 to 0422 

L3k 

2008 

278-293 

4 Oct- 19 Oct 

15 

91 

1283 to 0145 

L2d 

2008 

330-352 

25 Nov- 17 Dec 

22 

91 

0096 to 0423 

L2e 

2009 

068-101 

9 Mar- 11 Apr 

33 

91 

1286 to 0424 

L2f 

2009 

273-284 

30 Sep- 11 Oct 

11 

91 

1280 to 0084 


1 There are 119 tracks in the 8-day orbit and 1354 tracks in the 91 -day orbit. Cycle 
numbers are included for the 8-day repeat periods. 


Table A.2: Significant ICESat Parameters and Events by Campaign 


Cam- 

paign 

Year 

Day 

of 

year 

S/C 

orient- 

ation 1 

Start 

Beta’ 

Angle 

o 

End 

Beta’ 

Angle 

o 

Start 

Laser 

Infrared 

Energy 

(mJ) 

End 

Laser 

Infrared 

Energy 

(mJ) 

Mean 
footprint 
major 
axis (m) 

Day of year - 
comments 

- 

2003 

013 

- 

- 

- 

- 

- 

- 

013 - launch 

Lla 

2003 

051-088 

-Y/+X 

-45 

-32 

72 

51 

149 

080 - yaw flip 
085 - safe hold, 
adjust temperature 

L2a 

2003 

268-277/ 

277-322 

+Y 

51 

69 

80 

55 

100 

277 - orbit change 
286 - laser 
temperature anomaly 
287, 302 - adjust 
temperature 
311 - GPS solar 
flare anomaly 

L2b 

2004 

048-081 

+Y 

54 

40 

57 

33 

90 


L2c 

2004 

139-173 

-X 

13 

-4 

33 

5 

88 

142-147 - adjust 
temperature 

L3a 

2004 

277-313 

-Y 

-48 

-58 

67 

62 

56 

293 - adjust 
temperature 

L3b 

2005 

048-083 

-Y 

-56 

-45 

68 

54 

80 

054 - suspected 
amplifier bar drop, 
begin footprint 
anomaly 2 
068 - suspected 
amplifier bar drop 

L3c 

2005 

140-174 

+X 

-20 

-4 

49 

44 

55 


L3d 

2005 

294-328 

+Y 

51 

63 

43 

39 

52 


L3e 

2006 

053-087 

+Y 

62 

48 

38 

30 

52 


L3f 

2006 

144-177 

-X 

20 

4 

30 

30 

51 

149 - Energy jump 
up 2m J 

L3g 

2006 

298-331 

-Y 

-44 

-54 

30 

24 

53 

3 10 -begin ITRF 
2005 

L3h 

2007 

071-104 

-Y 

-60 

-47 

24 

21 

56 


L3i 

2007 

275-309 

+Y 

32 

46 

22 

20 

57 


L3j 

2008 

048-081 

+Y 

74 

62 

20 

16 

59 


L3k 

2008 

278-293 

+X 

-28 

-32 

18 

12 

52 

289 - Energy drop 4 
mJ 

L2d 

2008 

330-352 

-Y 

-45 

-53 

8 

4 

- 

343-344 - adjust 
temperature, energy 
up 5 mJ 

L2e 

2009 

068-101 

-Y 

-71 

-59 

6 

2 

- 

094-095 - adjust 
temperature 

L2f 

2009 

273-284 

-X 

20 

25 

4 

2 

- 


- 

2010 

242 

- 

- 

- 

- 

- 

- 

242 - reentry 


1 The spacecraft is said to be in “Sailboat” mode for ±Y orientations and in “Airplane” 
mode for ±X orientations, where the direction indicates the solar panel orientation with 
respect to the spacecraft velocity using the GLAS coordinate frame. 


2 The footprint diameter during L3b changed from a mean of 54 m (day of year 048-053) 
to 84 m (055-068). The reason for the larger footprint size during the latter part of the 
campaign is unknown, although a suspected amplifier bar dropout occurs near the event. 
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Appendix B Off-Nadir Atmospheric Delay Corrections Addendum 


The atmospheric delay calibration ATBD discussed the effects of off-nadir pointing on 
the atmospheric delay correction and concluded that for pointing within 15° of nadir, a 
simple l/sin(elevation angle) formulation would be provide a sub-millimeter accuracy 
corrections. Since it now seems likely that ICESat might point as far as 35o off-nadir, we 
have revaluated this approximation for large off nadir angles (up to 35°). The geometry 
for relating off-nadir angle, 0, to zenith angle, z=90-elevation angle (s), is shown in 
Figure 1. 


GLAS 



Figure 1: Geometry relating nadir angle to elevation angle. R is the radius of the Earth 
(~6378 km); h is the altitude of the satellite (~600km); 0 is the nadir angle; and z is the 
zenith angle, which equals 90-s, where e is the elevation angle. For q=35°, z=38.868° and 
e=5 1.132°. 

The comparison here we used the Niell [1996] hydrostatic mapping function that relates 
the atmospheric delay in the zenith direction to the delay at a specified (in- vacuum) 
elevation angle. The form of the mapping function is a continued fraction in sin(e). With 
the coefficients appropriate for polar regions the Niell mapping function takes the form: 


A L = m{£)lS.L z 


( 1 ) 


m(e ) = 


1/(1 + a/(l + bj\ + c))) 
sin(£) + a/(sin(£) + b/( sin(£) + c))) 


where AL is the atmospheric delay at elevation angle, s; AL Z is the delay in the zenith 
direction. For Polar Regions, under average conditions, a= 1.2046x1 O' 3 , b=2.90249xl0' 3 , 
and c=64.258xl0' 3 . 

Figure 2, shows the values of AL, for AL z =2.3m as a function of off-nadir angle for m(s) 
given in equation (1) and m(s) given simply by l/sin(e). At this scale, the differences are 
difficult to see. In Figure 3, we show the difference in units of mm. At the largest off- 
nadir angle, the difference is <2.5 mm and well within the atmospheric delay model error 
budget, 



Figure 2: Atmospheric delay as function of off-nadir angle under nominal conditions 
given in the text. Units are meters. 



Nadir Angle (deg) 

Figure 3: Difference between the Niell mapping function and the cosecant mapping 
function under nominal conditions. Difference units are millimeters. 

Bending effects 

The effects of bending were approximately evaluated in the ATBD for off-nadir angles 
up to 15o. We have more carefully considered these effects here because the off-nadir 
angles could be as large as 35°. To evaluate the effect we ray-traced through a standard, 
spherically symmetric, atmospheric model keeping careful track of the bending angles 
and the deviation between the vacuum and refracted paths. The ray tracing was 
performed from the ground to the satellite (at 600 km altitude) since this approach tends 
to more numerically stable than ray tracing from vacuum into the Earth’s atmosphere. 

The ray tracing started at a series of elevations ranging between 90 and 50 degrees. From 
the ray-trace, the nadir angle at the satellite and the angle subtended by the arc between 
initial starting point and the position of ray when it reached 600 km were computed. In 
addition, we also integrated the atmospheric delay and the bending angle as checks on the 
ray-trace. The differences between the atmospheric delays computed from the ray-trace 
and those given by the simply cosecant law were indistinguishable from those shown in 
Figure 3 above. The bending angle matched the values given be the Astronomical 
Almanac [1999] formula in the ATBD to within one milli-degree. 

Figure 4 show the arc distance to the footprint from the sub-satellite point computed from 
the ray tracing and from simple in vacuum geometry. Figure 5 shows the difference. In 
the worst case, the difference in foot print location is less than 5 meters and we conclude 
that atmospheric bending effects on foot print location can be ignored even for the largest 


off nadir pointing angles. Intuitively these results make sense when it is realized that 
most of the bending occurs in the lowest 10km of the atmosphere, and for a bending 
angle of 0.014 degrees from 10 km altitude, the foot print displacement would be 3.2 m 
for a ray at 50° deg elevation angle. 

Figure 4 also shows that for a 94° inclination orbit (sub satellite point ~440 km from the 
pole, that a nadir angle of ~35° will be needed to range to the pole. 

Conclusions 

Based on these calculations, we recommend that a simple l/sin(e) where e is given by 
cos' 1 (sin(0)R g /R s ) and R g is the radius at the footprint and R s is the radius to ICESat. To 
<5 meters, the bending of the ray by the atmosphere can be neglected in geo-locating the 
footprint. 



Nadir Angle (deg) 


Figure 4: Distance between the foot print location and the sub-satellite point as a 
function of nadir angle computed in vacuum and by ray tracing through a standard 
atmospheric delay model. 



Figure 5: Difference between the two curves in Figure 4 shown in meters. 
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